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Abstract 



In this paper we have tested several general numerical methods in solv- 
ing the quantum billiards, such as the boundary integral method (BIM) and 
the plane wave decomposition method (PWDM). We performed extensive 
numerical investigations of these two methods in a variety of quantum bil- 
liards: integrable systens (circles, rectangles, and segments of circular annu- 
lus), Kolmogorov-Armold-Moser (KAM) systems (Robnik billiards), and fully 
chaotic systems (ergodic, such as Bunimovich stadium, Sinai billiard and car- 
diod billiard). We have analyzed the scaling of the average absolute value of 
the systematic error AE of the eigenenergy in units of the mean level spacing 
with the density of discretization b (which is number of numerical nodes on 
the boundary within one de Broglie wavelength) and its relationship with the 
geometry and the classical dynamics. In contradistinction to the BIM, we find 
that in the PWDM the classical chaos is definitely relevant for the numerical 
accuracy at a fixed density of discretization b. We present evidence that it is 
not only the ergodicity that matters, but also the Lyapunov exponents and 
Kolmogorov entropy. We believe that this phenomenon is one manifestation 
of quantum chaos. 

PACS numbers: 05.45.+b, 02.70.Rw, 03.65.Ge 
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I. INTRODUCTION 



It is quite embarrassing to realize that in an attempt to numerically solve the Helmholtz 
equation 

V^(r) + k 2 ^{r) = 0, (1) 

satisfied by the scalar solution ip(r) with eigenenergy E = k 2 inside a connected plane do- 
main B with the Dirichlet boundary condition ^(r) = on the boundary dB, one can face 
enormous difficulties in cases of " problematic " geometries such as various nonconvex shapes. 
This is precisely the problem of solving and describing the quantum billiard B as a Hamil- 
tonian dynamical system, which is thus just the two-dimensional Schrodinger problem for a 
free point particle moving inside the enclosure dB, described by the wave function ip(r) with 
the eigenenergy E = k 2 . The corresponding classical problem is the classical dynamics of a 
freely moving point particle obeying the law of specular reflection upon hitting the boundary 
dB. Quantum billiards and their correspondence to their classical counterparts, especially 
in the semiclassical level, are important model systems in studies of quantum chaos |1] -Q. 
There are several general methods for a numerical solution of Eq.(|l|) such as the boundary 
integral method (BIM) (see e.g. @-§|), and the plane wave decomposition method (PWDM), 
introduced and advocated by Heller 0, whose analysis, especially in the light of the rele- 
vance of classical chaos, is the subject of our present paper. Another quite general method 
is the conformal mapping diagonalization technique introduced by Robnik || and further 
developed by Berry and Robnik ||, Prosen and Robnik [ID|, and Bohigas et al. which 



in principle should work for any shape, whereas, in practice it is used for shapes for which 
the conformal mapping onto the unit disk (or some other integrable geometries admitting 
a simple basis for the representation) is sufficiently simple (possiblly also analytic). These 
methods can face quite similar problems in cases of almost intractable geometries, but they 
are to some extent complementary. For example, the conformal mapping diagonalization 
technique can provide a complete set of all eigenenergies up to some maximal value beyond 
which the calculations cannot be performed due to the lack of computer storage (RAM), 
which means that we cannot reach very high-lying eigenstates. (Our present record 
is about 35 000 for the size of the banded matrix that we diagonalize in double precision, 
yielding at least 12 000 good levels with accuracy of at least 1% of the mean level spacing.) 

However, using the PWDM it is possible to go higher in energy by orders of magnitude 
but then only a few selected states can be calculated with many intermediate states in 
the spectral stretch missing. Therefore, the geometry of some interesting and representative 
high-lying states can be analyzed, but the sample is typically not sufficiently complete (there 
are many states missing) to perform statistical analysis. See, e.g., our recent papers on this 
topic | ]T2"| , |T3"| ] . The reasons for a failure of one of these methods can be quite different. For 
example, in the BIM the main difficulty stems from the existence of " exterior chords" in 
nonconvex geometries in its standard formulation (see Sec. Ill), but the trouble might be 
overcome by an appropriate reformulation of the method adapted to the correct semiclassical 
behaviour. We will discuss this in Section III where we also show that classical chaos 
is completely irrelevant for the BIM. On the contrary, in the PWDM we find that the 
classical chaos is relevant for numerical accuracy especially in the semiclassical limit of the 
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sufficiently small effective Planck constant h e ff reached at sufficiently high eigenenergies. 
This demonstration and its qualitative explanation is the main subject of our present paper. 
To give a specific example we should mention isospectral billiards discovered and proved 
by Gordon et al. ||14|| , which have been investigated experimentally by Sridhar and Kudrolli 



|T5| and it is also our experience that the BIM fails in this case (namely, due to strong 
nonconvexities), whereas the PWDM at b = 12 yields the accuracy of eigenenergies of 
within a few percent of the mean level spacing, except for some very special eigenmodes 
for which surprisingly we find agreement within double precision (16 digits) and which are 
characterized by the fact that these eigenvalues agree with the analytic solutions for the 
triangles within single precision (eight digits). So the fact that in this and similar cases the 
experimental precision (for some levels) exceeds the best possible numerical precision even 
when using the best available methods is embarrassing for a theoretician, but also motivation 
for further work. 

The paper is organized as follows. In Sec. II, we focus our attention on the PWDM. In 
Sec. Ill, we first point out some serious flaws in the derivation of the BIM in the literature 
and show how the final formula (which nevertheless was correct) should be derived in a 
regularized way and then discuss the numerical results of the BIM and the relevance with 
classical dynamicss and geometry, etc. In Sec. IV we give a discussion and conclusions. 



II. THE PLANE WAVE DECOMPOSITION METHOD OF HELLER 

A. The numerical procedure of the PWDM 

In this section we present our general exposition of PWDM following |12| . To solve the 
Schrodinger equation ([!]) for ^>(r) with the Dirichlet boundary condition ip(r) =0 on dB 
we use the ansatz of the superposition of plane waves (originally due to Heller J7[) 

TV 

^( r ) = H a i c os(k xj x + k yj y + fa), (2) 

3=1 

where k x j = kcosdj, k y j = ksin9j, k 2 = E, and we use the notation r = (x,y). N 
is the number of plane waves and <pj are random phases, drawn from the interval [0,27r), 
assuming a uniform distribution, and 8j = 2jn/N determining the direction angles of the 
wave vectors chosen equidistantly. The ansatz (Q) solves the Schrodinger equation (|l|) in 
the interior of the billiard region B, so that we have only to satisfy the Dirichlet boundary 
condition. Taking the random phases, as we discovered, is equivalent to spreading the origins 
of plane waves all over the billiard region, and at the same time this results in reducing the 
CPU time by almost a factor of 10. For a given k we set the wave function equal to zero at a 
finite number M of boundary points (primary nodes) and equal to 1 at an arbitrarily chosen 
interior point. Of course, M > N. This gives an inhomogeneous set of equations that can 
be solved by matrix inversion. Usually the matrix is very singular and thus the singular 
value decomposition method has been invoked |T|,|T7f . After obtaining the coefficients aj we 
calculate the wave functions at other boundary points (secondary nodes). We always have 
three secondary nodes between a pair of primary nodes. The experience shows that a further 
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increase of the number of secondary nodes does not enhance the accuracy. The sum of the 
squares of the wave function at all the secondary nodes (Heller called this sum " tension " ) 
would be ideally zero if k 2 is an eigenvalue and if Eq.(^) is the corresponding exact solution 
of Eq.(JT]). In practice it is a positive number. Therefore, the eigenvalue problem now is 
to find the minimum of the tension. In our numerical procedure we have looked for the 
zeros of the first derivative of the tension; namely, the derivative is available analytically or 
explicitly from Eq.(^) once the amplitudes aj have been found. In this paper we make the 
choice M = N since it proved to be sufficient for calculating the lowest 100 states whose 
accuracy we analyze. (For the high- lying states studied in our paper ]T2|, we have used 
M = 5N/3.) It must be pointed out that the wave functions obtained in this way are not 
(yet) normalized due to the arbitrary choice of the interior point where the value of the wave 
function has been arbitrarily set equal to unity. We therefore explicitly normalize these wave 
functions before embarking on the analysis of their properties. 

The accuracy of this method of course depends on the number of plane waves (N) and 
on the number of the primary nodes (M) and we have a considerable freedom in choosing 
N and M > N. In order to reach a sufficient accuracy the experience shows that we should 
take at least N = 3C/\d e Brogiie and M = N, where C is the perimeter of the billiard and 
Ade Brogiie is the de Brogiie wavelength equal to 2n/k. With this choice in the present context 
and for the lowest 100 states we reach the double precision accuracy (16 digits) for all levels 
of integrable systems such as the rectangular billiard (where the eigenenergies can be given 
trivially analytically) and the circular billiard, but also for the Robnik billiard B\ for small 
A < 0.1. Introducing the density of discretization b defined as the number of numerical 
nodes per one de Brogiie wavelength on the boundary we thus write the number of plane 
waves N = bC/\ de Brogiie = b2n£/k. 

The main problem of investigation in this paper is to study the dependence of the sys- 
tematical numerical error AE (i.e., the error due to the finite discretization) on the density 
of discretization b and the dependence of AE on the geometry (billiard shape parameter) 
at fixed b. In order to perform a systematic analysis the errors should be measured in some 
natural units and in our case this is of course just the mean level spacing, which, accord- 
ing to the leading term of the Weyl formula, is equal to Air/ A, where A is the area of the 
billiard B. From now on we shall always assume that AE of a particular energy level is in 
fact measured in such natural units. Of course one immediately realizes that the error AE 
fluctuates wildly from state to state (see Figs. 5 and 7) so that generally nothing can be 
predicted about it individually. Therefore, the approach must be a statistical one and so we 
typically take an average of the errors AE over a suitable ensemble of states. Specifically, 
in all cases of this paper we have taken the average of the absolute values of AE over the 
lowest 100 states (of a given symmetry class) and denoted it by (\AE\). It is important and 
should be mentioned that we have also checked the stationarity of such an average value over 
consecutive spectral stretches of 100 states each, so that our procedure does make sense. In 
addition, we have also investigated the standard deviation <j\ae\ — ((\AE\ — (\ AE])) 2 ) 1 / 2 , 
which always has the same order of magnitude as the average value. 

It turns out that the accuracy of energy levels depends nontrivially on b, unlike in the 
BIM, where we find always a power law (see Sec. Ill), namely, it typically shows broken 
power law. By this we mean that (|AJ5|) obeys a power law 



4 



(\AE\) = Ab~ a , 



(3) 



with very large a for sufficiently small 6, b < b c , whereas for larger b > b c it obeys a rather 
flat power law with very small positive a (close to zero). Therefore, in contradistinction 
to the BIM, it is difficult to explore the general dependence of (|A£J|) on b if there is any 
such universality at all. However, in order to investigate the dependence of the accuracy on 
geometry and the implied dynamical properties of billiards, we have decided to fix the value 
of b and have chosen b — 12, and then we look at the dependence of (|A.E|) on the shape 
parameters of three one-parameter billiards, namely, the Robnik billiard, the Bunimovich 
stadium, and the Sinai billiard. 

Finally, we would like to discuss how to estimate the error. As usual, to speak of an error 
we need to have a standard value. The question is how to get the standard values in different 
quantum billiards. As for the integrable billiard, such as a rectangle, we know analytically 
the exact values. For a circular billiard, they are the zeros of the Bessel function, which can 
be calculated very precisely. However, for other billiards, in particular the chaotic billiards, 
there are no true accurate values available (otherwise we do not need the numerical methods 
anymore). In fact, both the PWDM and the BIM can be self-tested for their accuracy. On 
the one hand, in both cases, the numerical value at very large b can be regarded as the " true " 
value. On the other hand, in the PWDM one may change the position of the interior point 
and compare the two lists of eigenenergies obtained. Moreover, since we have also some other 
special methods invented for the billiard of a specific geometry, such as the diagonalization 
method for the Robnik billiard and the scattering approach for the Sinai billiard, in these 
techniques the accuracy is well controlled and we may obtain more accurate results than the 
BIM and PWDM; thus we can use the eigenenergy list from them as the standard value. 

In our studies in this paper, the " standard value " of the Robnik billiard are provided by 
Prosen []T5|| by using the diagonalization method with a very large dimension of the matrix 
and thus the lowest 1000 eigenvalues are guaranted with an accuracy of at least 10~ 12 in 
units of the mean level spacing for the large shape parameter A. The eigenvalues of the 
Sinai billiard were provided by Schanz and Smilansky |19| by using their scattering method. 
The accuracy is about 10~ 7 of the mean level spacing, which is already high enough for 
our purpose. For the billiards whose eignevalues are not available from other methods, we 
always take the eigenvalues at very large density of discretization b (say 30) as the true 
value. 



B. Relevance of chaos with the numerical accuracy 

As is well known, in the classically integrable quantum Hamiltonian systems in the 
semiclassical limit (of sufficiently small h) the eigenfunction can locally be described by a 
finite superposition of plane waves with the same wave number; in the case of plane billiards 
it is k = \/E. If the quantum system has ergodic classical dynamics then in the semiclassical 
limit locally the wave function can be represented as a superposition of infinitely many plane 
waves with the same k and with the wave vectors being isotropically distributed on the 



circle of radius k |20]. Moreover, the ergodicity suggests to assume random phases for the 



ensemble of plane waves, which implies that to the lowest approximation the wave function 
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is a Gaussian random function. While this is a good starting approximation, originally due 



to Berry p0[ and recently verified by Aurich and Steiner |21] and also by Li and Robnik |T2 
the phases are actually not random but correlated in a subtle way dictated by the classical 
dynamics, especially along the short and the least unstable periodic orbits, which is the 
origin of the scar phenomenon p|j2^-p5|. Thus we can qualitatively very well understand 
that the PWDM should work well or even brilliantly in cases of classically integrable billiards 
whereas in the ergodic systems we expect a severe degradation of the accuracy (at fixed 
b) simply because the finite number of plane waves cannot capture the correct (infinite) 
superposition of plane waves everywhere in the interior of the billiard. If the system is a 
generic system of a mixed type with regular and irregular regions coexisting in the classical 
phase space, a scenario described by the Kolmogorov-Arnold-Moser (KAM) theory, then 
the degradation of accuracy (at fixed b) with increasing fractional measure of the chaotic 
component (denoted by P2) is certainly expected. However, p 2 is not the only parameter 
that controls the accuracy (at fixed b) since, as we shall see, the dynamical properties such 
as the diffusion time, Lyapunov exponent, and Kolmogorov entropy also play a role. It is 
the aim of the present paper to numerically explore this type of behavior in three different 
billiard systems. 

The first billiard system is defined as the quadratic (complex) conformal map w = z+Xz 2 
from the unit disk \z\ < 1 from the z plane onto the w = (x, y) complex plane. The system 



has been introduced by Robnik [26[ and further studied by Hayli et al. ||27|| , Frisk [28 
and Bruus and Stone |29[ for various parameter values A. Since the billiard (usually called 
the Robnik billiard by the people in the community) has an analytic boundary it goes 
continuously from the integrable case (circle, A = 0) through a KAM-like regime of small 
A < 1/4 with mixed classical dynamics and becomes nonconvex at A = 1/4 (the bounce map 
becomes discontinuous), where the Lazutkin caustics (invariant tori) are destroyed giving 
way to ergodicity. It was shown by Robnik that the classical dynamics at these values 
of A is predominantly chaotic (almost ergodic), although Hayli et al. |S7| have shown that 
there are still some stable periodic orbits surrounded by very tiny stability islands up to 
A = 0.2791. At larger A we have reason and numerical evidence [ISO] to expect that the 



dynamics can be ergodic. Recently, it has been proven rigorously by Markarian |3T[ that 
for A = 1/2 (a cardioid billiard) the system is indeed ergodic, mixing, and K. This was a 
further motivation to study the cardioid billiard classically, semiclassically, and quantanlly 
by several groups, e.g., Backer et al. [^2|, and Bruus and Whelan |j33| . The billiard shape 



for A = 0.4 is shown (the upper half) in Fig. 1(a). Since all states are either even or odd we 
can take into account these symmetry properties explicitly. In fact, we want to specialize 
to the odd eigenstates only. Therefore, in order to a priori satisfy the Dirichlet boundary 
condition on the abscissa of Fig. la we specialize the general ansatz (El) to the form 



N 



^( r ) = a j cos(k xj x + sm(k yj y), (4) 
j=i 

where all the quantities are precisely as in Eq. (EJ) except that the N discretization (primary) 
nodes are equidistantly located only along the half of the full billiard boundary, so that b is 
exactly the same as in using the ansatz (EJ) for the full billiard. 

In Fig. 2 we show the results for this billiard, namely we plot (\AE\) (in logarithmic 
units) versus A at fixed b = 12. Close to integrability (A < 0.1) we reach the accuracy within 
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14-15 digits, which is almost the double precision on our machine (16 digits), in which all 
our calculations have been performed. As the value of A increases we observe a dramatic 
deterioration of the accuracy where (|AJ5|) increases by many orders of magnitude, namely, 
by almost 13 decades, leveling off at (\AE\) approximately equal to 10~ 2 , which means that 
we have now the accuracy of only a few percent of the mean level spacing. This dramatic 
but quite smooth increase of (|Ai£|) is certainly related to the emergence of classical chaos 
with increasing A, but definitely is not controlled merely by p2 because p2 reaches the value 
of 1 (almost ergodicity) already at A = 1/4 []T0| , p6|| , whereas (|A.E|) still varies considerably 



in the region A > 1/4. Thus it is obvious that in the semiclassical picture also other 
classical dynamical properties (measures of the " hardness " of chaos) play an important 
role. Although we do not have a quantitative theory yet one should observe that according 
to Robnik |26[] the Lyapunov exponent and Kolmogorov entropy (h) vary also quite smoothly 
with A, suggesting a speculation that there might be a relation between (|A£?|) and h. 

Another demonstration of the effectivity of the PWDM and its accuracy is displayed 
in Table I where we show the numerical value of the scalar product of two consecutive 
normalized eigenstates, namely, the ground state and the first excited state, denoted by 
O12, which ideally should be zero. We see here too that the accuracy decreases (by orders 
of magnitude) sharply but smoothly with increasing shape parameter A. 

It is then interesting to similarly analyze an ergodic system such as the stadium of 
Bunimovich shown in Fig. 1(b), where the shape parameter is a/R and we have looked at 
the results for < a/R < 10. In fact, for our purposes we have chosen and fixed R = 1 in 
all cases. Since this billiard is known to be rigorously ergodic (and mixing and K) for any 
a > in this case p2 is exactly 1 and constant. We have calculated the lowest 100 energy 
levels of the odd-odd symmetry class. Therefore, in this case the general ansatz (||) can be 
specialized as 

N 

VK r ) = a i sin(fc a;i a;) sin(k yj y). (5) 

3=1 

Here again the discretization (primary) nodes are only on the outer boundary of the stadium 
with discretization density b = 12. From our plot in Fig. 3 we see that in the integrable 
case of the circle (a = 0) we again reach the accuracy within at least 14 digits, but this 
brilliant accuracy at fixed b = 12 deteriorates almost discontinuously upon increasing a and 
then (|A.E|) still increases by about two orders of magnitude when a goes from 0.1 to 10. 
It appears to us that classical chaos is definitely relevant for the accuracy of the method 
which might and should be explained by an appropriate theory in the semiclassical level. As 
an observation we should mention that the Kolmogorov entropy increases sharply with a/R 
when a/R goes from to about 1, where it reaches the maximum, and then decreases slowly 
34"], whereas our (\AE\) increases monotonically. Thus if there is a relationship between 



(|Ai<7|) and Kolmogorov entropy it certainly is not a simple one. 

We have tested also another system with hard chaos, namely, the Sinai billiard sketched 
in Fig. 1(c) (desymmetrized). The system is known to be ergodic, mixing, and K. In 
calculating the 100 lowest energy levels of the desymmetrized Sinai billiard we used the 
same specialized ansatz as in Eq. (|j), thereby taking into account explicitly the Dirichlet 
boundary condition on the abscissa y = 0. In this case b is the density of discretization of 
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the equidistant nodes along the rest of the perimeter. Similarly as in the case of the stadium 
we easily reach the double precision of 16 digits in the limiting integrable case of zero radius 
R = 0, but this accuracy is almost instantly lost by increasing R, as seen in Fig. 4. The 

value of (\AE\) levels off at about 10" 4 10~ 2 for all R between 0.025 and 0.45. 

As a final technical point we comment on the stationarity of (|A.E|) as a function of 
energy, which has been confirmed for the Robnik billiard at A = 0.27, where the average 
value over consecutive spectral stretches over 100 states has been found to be quite stable. 
Specifically, to illustrate this finding we plot in Fig. 5 the absolute values of the errors of 
the lowest 400 consecutive eigenstates where one can see that the average value over 100 
consecutive states is quite stable indeed. This is shown in Table II for four intevals of 100 
states each. We should emphasize again that the fluctuation is very big. We have calculated 
the standard deviation for these 400 levels; it is <j\ae\ — 7.40 x 10~ 7 , which has the same 
order of magnitude as the average value (see Table II). It is our numerical experience that 
for all cases that the standard deviations are always about the same order of average values. 



III. BOUNDARY INTEGRAL METHOD 

After discussing the PWDM, we now turn to another very important numerical method: 
the boundary integral method. This method is widely used not only in studying quantum 
chaos, but also in engineering [KB]. In this section we would give an extensive numerical 



investigation of its accuracy in relation to classical dynamics and geometrical properties. 
However, before we go into a detailed numerical and technical analysis, we would like to 
point out two serious flaws in the derivation of the BIM in the literature and show how the 
final formula, which nevertheless is correct, should be derived in a sound way. 



A. A regularized derivation of the BIM 

In order to clearly expose the difficulties and the errors in the derivation of the BIM 
offered in the literature, e.g. in Ref [||], we present our regularized derivation, by which 
we mean that we construct and use a Green function that automatically (by construction) 
satisfies the Dirichlet boundary condition (vanishes locally on the boundary dB of the billiard 
domain B), which is achieved by employing the method of images (see, e.g., the article of 
Balian and Bloch and the references therein). This will enable us to avoid committing 



two errors, which, however, luckily compensated for each other. First, in taking the normal 
derivatives on the two sides of Eq. (6) in Ref , on the right hand side we must use the value 
^(r) which is the interior solution inside B, rather than ^ip(r), which is the value exactly 
on the boundary, simply because in taking the derivatives we must evaluate the function at 
two infinitesimally separated points normal to the boundary. Second, this error of taking 
the unjustified factor 1/2 is then exactly compensated for by another error in arriving at 
Eq. (8) in Ref. 0, namely, by interchanging the integration along the boundary dB and the 
normal differentiation, because due to singularities on the boundary these two operations 
do not commute. 

Now, we offer our regularized derivation. We are searching for the solution ip(r) with 
eigenenergy E = k 2 obeying the Helmholtz equation (|1|), with the Dirichlet boundary con- 
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dition ^(r) = on the boundary r G dB. We will transform this Schrodinger equation for 
our quantum billiard B into an integral equation by means of the regularized Green function 
G(r, r'), which solves the defining equation 

V 2 r G(r, r') + k 2 G(r, r') = 5(r - r') - 5(v - r' R ), (6) 

where r and r' are in B U dB and r'^ is the mirror image of r' with respect to the tangent at 
the closest-lying point on the boundary (and thus if r' is sufficiently close to the boundary 
then r'^j is outside the billiard B) (see Fig. 6). The solution can easily be found in terms of 
the free propagator (the free-particle Green function on the full Euclidean plane) 

G (v,r') = ~iH^(k\r-r'\), (7) 

where Hq is the zeroth order Hankel function of the first kind |37| , namely, 

G(r,r') = Go(r,r')-Go(r,r' R ), (8) 

such that now G(r, r') is zero by construction for any r' on the boundary dB, in contradis- 
tinction to the Green function defined and used in Eq. (5) in ||. Multiplication of Eq.(^|) 
by ip(r) and the Helmholtz equation ([!]) by G(r, r'), subtraction, integration over the area 
inside B, and using Green's theorem yields 

j> ds [V>(r)n ■ V r G(r, r) - G(r, r')n • V r ^(r)] = ^(r'), (9) 

where s is the arclength on the boundary dB oriented counterclockwise, n is the unit normal 
vector to dB at r oriented outward, and this equation is now valid for all r' inside and on 
the boundary of B. Since in this equation everything is regular, we can take the normal 
partial derivatives on both sides. Following the usual notation in || we define the normal 
derivative of ip at the point s as 

u(s) = n- V r V>(r(s)) (10) 
and thus using the boundary condition ip(r) = we arrive at 

u(s) = -2 £ ds'u(s')n ■ V r G (r, r'). (11) 

In this way we have correctly derived the main integral equation of the boundary integral 
method, which is correctly given as Eq. (8) in (where the two errors exactly compensate 
for each other), so that all the further steps in working out the geometry of Eq. (|11|) and 
the numerical discretization are exactly the same as in 0]. As shown in Fig. 6, we define 
the length of the chord between two points on the boundary r(s) and r(s') as 

p(s,s') = \r(s)- r (s')\ (12) 

and the angle 6(s', s) is the angle between the chord and the tangent to dB at s'. The chord 
is an oriented separation vector pointing from r(s) to r' = r(s'). Of course 9(s, s') ^ 9(s', s). 
Thus, following the notation of || we can write 
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(13) 



and we obtain finally 

u (s) = --ik j ds'u(s') sinfl (s, a^HP {kp{s, s')}. (14) 

In numerically solving this integral equation we have used precisely the same primitive 
discretization procedure as [||] , which turned out to be better than some other more sophis- 
ticated versions. So we simply divide the perimeter £ into N equally long segments and 
thus define 

s m = mC/N, p(si, s m ) = pi m , 9(si, s m ) = 9i m , 1 < I, m < N. (15) 

Therefore, numerically, we are searching for the zero of the determinant An(E) =det(M; m ), 
where M/ m are the matrix elements of the N x N matrix 

Mi m = 5i m + — — sin O lm H?\k Plm ), (16) 

where E = k 2 . Due to the asymmetry 6i m ^ 6 m i this matrix is a general complex non- 
Hermitian matrix. For the diagonal elements / = m, where Q\ m is either zero or tt, the 
proper limit of the second term on the right hand side in Eq. (|16|) must be taken and then 
it is equal to k(s)£/2iiN , where k(s) is the curvature of the boundary at s. 

One important aspect of this formalism is the semiclassical limiting form that has been 
extensively studied by Boasman ||. At this point we want to make the following comment. 
In the cases of nonconvex geometries we will have exterior chords connecting two points on 
the boundary such that they lie entirely, or at least partially, outside B. While formally 
the method and the procedure in such cases is perfectly correct, in reality it might be prob- 
lematic, which can be seen by considering the semiclassical limit. The formal leading order 
in the asymptotic expansion of the Hankel function H^' (Debye approximation) does not 
match the actually correct semiclassical leading approximation that would be spanned by the 
shortest classical orbit connecting the two points via at least one or many collision points in 
between. Therefore, we understand and expect that the method must meet some difficulties 
in cases of nonconvex geometry. This has been partially confirmed in our present work as 
we will show in Sec. IIIB, while the analytical work to reformulate the method including the 
multiple collision expansions is in progress |38| and is expected to deal satisfactorily with 



nonconvex geometries. 

In the Sec. IIIB we shall analyze the numerical accuracy of the BIM as a function of the 
density of discretization 

. 2ttN 



k£ 



(17) 



in a variety of quantum billiards with integrable, KAM-type, or ergodic classical dynamics, 
including such with nonconvex geometry. The main result is that there is always a power 
law so that the error of eigenenergy in units of the mean level spacing, after taking the 
average of the absolute value over a suitable ensemble of eigenstates, obeys (|Ai?|) = Ab~ a , 
but the exponent a (and the prefactor A) is nonuniversal. 
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B. Numerical results 



The numerical procedure we have used to solve the variety of quantum billiards is exactly 
as described above and therefore it is precisely the same as in |J. Our main task now is 
to analyze in detail the behavior of the BIM as a function of the density of discretization 
b, especially in relation to the geometrical properties of B (nonconvexities) and in relation 
to classical dynamics, whose chaotic behavior is expected to imply interesting methodologic 
and algorithmic manifestation of quantum chaos. 

Again, like in the PWDM, we have to measure the numerical error of the eigenenergies in 
units of the mean level spacing and perform some kind of averaging over a suitable ensemble 
of consecutive states. However, this will make sense only if such a local average of the error 
is stationary (constant) over a suitable energy interval. This condition has been confirmed 
to be satisfied in almost all cases that we checked. In the case of the circle billiard and in the 
case of the cardioid billiard this stationarity of the locally averaged error is shown in Figs. 
7(a-b), respectively, where we plot the data for about 95% of the lowest 1000 odd levels. The 
average value (|A.E|) for these two cases are 3.89 x 10~ 5 ( circle) and 1.99 x 10~ 3 (cardiod 
billiard), while the standard deviations <j\ae\ are 3.67 x 10 -5 (circle) and 2.91 x 10~ 3 (cardiod 
billiard), respectively. Again, like in the case of the PWDM, the standard deviation in the 
BIM is of the same order as the average value. 

We have established that we have always taken the average of the absolute value of the 
error (in units of the mean level spacing) over a suitable ensemble of eigenstates, for which 
we have chosen the lowest 100 eigenstates in all cases. Actually, strictly speaking, we have 
taken about 90-95 levels from the ensemble of the lowest 100 odd levels: since in using 
the BIM we always miss some levels, the quota of missing levels typically is a few percent, 
depending on the step size. In our case the step size is 1/20 of the mean level spacing and 
the fraction of missing levels varied from about 10 percent in integrable billiard with Poisson 
statistics to about 5 % in ergodic cases with GOE statistics. 

In Figs. 8(a-f) we show the error (|A£?|) versus b for six billiard shapes. In Fig. 8(a) 
we have the full circle billiard. In Fig. 8(b) we have 270° (= 37r/2) segment of the annulus 
billiard with inner radius R\ = 0.45 and outer radius R2 = 0.5. These data were based 
on the very careful work of Thomas Hesse | 39| , who has kindly communicated to us his 
unpublished results and the analysis, which we have independently checked and confirmed. 
In Figs. 8(c) and (d) we have the Robnik billiard with shape parameters A = 0.15 and 
A = 1/2, respectively. In Fig. 8(e) we show the results for the 1/4 Bunimovich stadium 
with the size 2x2 for the central square. In Fig. 8(f) we have the 1/4 Sinai billiard with 
dimensions 2x2 for the square and circular radius R = 1/2. The best-fitting power-law 
curve is described by Eq.|3| and is seen to provide a very significant fit in all six cases. 

As for the Robnik billiard, one should be reminded that at A = we have integrable 
classical dynamics in the circle billiard and at A = 0.15 we have KAM-type dynamics 
with islands of stability fL6 |; it should be emphasized that at A = 1/4 we have a zero 
curvature point at z = — 1 and for all A > 1/4 the shape is nonconvex, whilst at A = 1/2 
we have ergodicity and also nonconvex geometry. Technically, in all cases at various A we 
have calculated all states by applying the BIM, but then, for technical reasons compared 
only the odd states with their exact value, which are supplied by the conformal mapping 



diagonalization technique [8,10 
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In Figs. 8(a-f) we thus observe that there is no clear relationship between the value of 
a and the degree of classical chaos in all the various billiards. However, it is an interesting 
" experimental " result that the power law (0) seems to be universally valid, with nonuniversal 
numerical value of a and A. 

Having established the validity of the power law (|3|) it is now most interesting and also 
immensely CPU time consuming (it took almost one month of CPU time on a Convex C3860 
to produce Figs. 9 and 10) to look at the variation of a with the billiard shape parameter 
A, which is shown in Fig. 9. There is a flat region of almost constant a within < A < 1/4: 
It fluctuates slightly around 3.5. At A > 1/4 the nonconvexities of the boundary appear; 
unlike naive expectation, a now even increases up to a value of slightly larger than 5 reached 
at A ~ 0.35 and then starts to decrease rapidly down to the value of about 2 at A = 1/2. 
Therefore, there is no clear correlation with either the nonconvexities or the classical chaos. 

In addition to a in Eq.@ we would also like to know the value of the constant A (the 
pre-factor) in each case. This is given in Fig. 10 by fixing b = 12 and plotting the mean 
absolute value of the error (averaged over the lowest 100 odd states) versus A. Here we 
see that the mean error (|A.E7|) is almost constant up to A < 0.35 and is equal to about 
6 x 10~ 6 . At A > 0.35 we observe the rapid increasing of (|A.E|}. Here again we cannot 
draw a clear conclusion but only note that the error starts to increase rapidly in the regime 
where classical chaos becomes " hard" (the Kolmogorov entropy increases steeply ||26||). 

Most of our results are summarized in Table III for three classes of billiard systems with 
different type of classical dynamics, namely, integrable, KAM-type, and ergodic systems. 
We show the calculated values of a and also the average absolute value of the error (\AE\) 
with fixed value of b = 12. The table clearly demonstrates that the power law (|3|) for 
the BIM is universal, but not the exponent a and the prefactor A. It also demonstrates 
that classical dynamics has little effect on a, whereas the nonconvexities of the boundary 
might be more important. There is no theory for a so far except for the circle billiard, for 
which Boasman || predicts a = 3, which might be marginally compatible with our data. 
The bizzare behavior of a in various dynamical regimes reminds us of the difficulties in 
theoretical predictions of, e.g., classical correlation functions p0| . In the table we include 
also the results of the integrable case of the rectangular equilateral triangle (half of the unit 
square) where a = 3.28 is quite large and the ergodic case of the 1/4 Heller's stadium (2x2 
square plus two semicircles with a unit radius) in which case a = 3.0 is also quite large. 

When thinking about improving the efficiency and the accuracy of the BIM we have 
also tried a more sophisticated version of the BIM, where we have explicitly used a Gaus- 
sian integration on the boundary when discretizing our main equation ([11]). However, this 
experience has been negative after many careful checks in various billiards and therefore 
we decided to resort to the primitive discretization of the BIM, which is exactly the same 
approach as in ||. 



IV. DISCUSSION AND CONCLUSIONS 



We believe that our present paper presents quite firm numerical (phenomenological) ev- 
idence for the relevance of classical chaos for the effectiveness of the PWDM as a quantal 
numerical method to solve a quantum billiard, which is manifested especially in the semi- 
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classical limit and might and should be explained in terms of an appropriate semiclassical 
theory. Qualitatively, the reasons for this phenomenon are explained before. The parameter 
P2, the fractional volume of the chaotic component (s), definitely plays an important role, 
but is not the only aspect of classical chaos controlling the behavior of the error (\AE\) at 
fixed discretization density b. Namely, even in rigorously ergodic systems where p2 = 1 the 
error (\AE\) might be controlled by the slow diffusion in the classical phase space (diffusive 
ergodic regime, soft chaos). If the classical diffusion time is much longer than the break time 
hreak (tbreak = fi/D, where D is the mean energy level spacing) then the quantal states 
will be strongly localized in spite of the formal ergodicity (for a demonstration see Ref. ITS 



for the Robnik billiard and Refs. [|l],fEJ for the stadium billiard) and therefore they mimic 
a certain amount of regularity, enabling a better accuracy of the PWDM, i.e., (\AE\) is 
smaller than for completely extended chaotic high-lying eigenstates where, according to our 
experience, somehow (\AE\) typically saturates at about a few percent of the mean level 
spacing, even if we drastically increase b beyond any reasonable limits. Indeed, as can be 
seen by a comparison of Figs. 2-4, in the case of the stadium this saturation value of (| A£7|)is 
about 10~ 4 , which is almost two orders of magnitude smaller than in the Sinai billiard (Fig. 
4) and the cardioid billiard (Fig. 2). We think that this is due to the strong localization 
of eigenstates in the stadium, which is very well known to display an unusual abundance of 
scars |7,23,25]. Thus our present work is a motivation for a semiclassical theory to explain 
this aspect of quantum chaos that exhibits some algorithmic properties of the PWDM in 
applying it to quantum billiards with a variety of classical dynamics. 

Further, we have investigated the behavior of the boundary integral method with respect 
to the density of discretization b as defined in Eq. (|T7| ) (6 is the number of numerical nodes 
per de Broglie wavelength along the boundary) since we expected some relevance of noncon- 
vexities and possibly of classical chaos. In all cases we discovered that there is a power-law 
behavior described in Eq. (Q). We wanted to verify whether there is any systematic effect 
of classical dynamics of quantum billiards on a and A. The answer is negative. On the 
other hand, we found that the role of nonconvex geometry of the boundary might be more 
important, although no final conclusion is possible at this point. The difficulties of the BIM 
might be expected as explained before, and the easiest way to see that is to consider the 
semiclassical limiting approximation of the BIM |[38|| . After explaining two systematic errors 
in the literature where the integral the BIM equation is derived and where luckily the two 
errors mutually compensate for each other exactly, we have given the correct (regularized) 
derivation and discussed the BIM formalism thus derived. We agree that even in nonconvex 
geometries it is formally right, but nevertheless practically might be less efficient, which is 
expected by considering the semiclassical limit mentioned above. Therefore, we suggest a 
generalization of the BIM by using a multiple reflection (collision) expansion in calculating 
the most appropriate Green function, which is another subject of our current investigation 
38| and is important not only for studies in quantum chaos but also in engineering prob- 



lems [[35[ since it would lead to a better efficiency of the BIM, namely, larger a. Most of our 
results are summarized in Table III, giving the evidence for the above conclusions. 

It remains an interesting and important theoretical problem to study the sensitivity of 
the eigenstates (eigenenergies and wave functions) on the boundary data of eigenfunctions, 
of which one aspect is also the dependence of the eigenstates on the billiard shape parameter. 
If such sensitivity correlates with classical chaotic dynamics and at the same time manifests 
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itself in the accuracy of the purely quantal numerical methods, then such a behavior would 
be one important manifestation of quantum chaos. This interesting line of thought in the 
search for another aspect of quantum chaos has been further developed in another work [43 



where we also present detailed studies of a level curvature distribution and other measures 
of the sensitivity of the eigenstates. 
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TABLE I. Test of the orthogonality of the eigenstates and the scalar product of two 
consecutive normalized wave functions O12, namely, the ground state and the first excited 
state, for the Robnik billiard at different shape parameters. The number in brackets repre- 
sents thepower of 10. 



A 


12 





2.0[-16] 


0.1 


-5.0[-15] 


0.2 


7.8[-10] 


0.3 


4.8 [-6] 


0.4 


5.5 [-4] 


0.5 


1.5 [-3] 



TABEL II. Stationarity test of (\AE\) for the Robnik billiard at A = 0.27 for the lowest 
400 odd eigenstates. 



Average stretch 


(|A£|> 


1-100 


1.54[-7] 


101-200 


2.21[-7] 


201-300 


2.77[-7] 


301-400 


2.03[-7] 


1-400 


2.13[-7] 
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TABLE III. Power-law exponent a and the average absolute value of the error (\AE\) 
with b = 12 for different billiards. For details of the KAM-type billiards see also Figs. 9 and 
10. 



Type 


Quantum billiard 


a 


(|A£|) b=12 


integrable 


circle (half) 
circle (full) 
rect angle-triangle 
segment annulus 


2.94 ± 0.17 
3.44 ± 0.18 
3.28 ± 0.29 
2.23 ± 0.24 


6. 74 [-5] 
5.97[-6] 
4.08[-5] 
1.77[-3] 


KAM 


Robmk (full) 

(0 < A < 1/4) 

Robnik (half) 

(0 < A < 1/4) 


« 3.4 
« 2.9 


« 5.0[-6] 
« 7.0[-5] 


ergo die 


stadium (1/4) 

cardioid (full) 

Sinai (1/4) 

Robnik (full) 
(0.3 < A < 1/2) 


3.00 ± 0.16 
2.10 ± 0.13 
2.47 ± 0.05 

see Fig. 9 


1.18[-4] 
3.04[-4] 
3.39[-3] 

see Fig. 10 
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FIGURES 



FIG. 1. Geometry of the boundary of the three desymmetrized billiards: (a) the Robnik 
billiard, (b) the Bunimovich stadium, (c) and the Sinai billiard. 

FIG. 2. Ensemble-averaged (over 100 lowest odd eigenstates) absolute error (measured in 
units of the mean level spacing) (\AE\) (in logarithmic units) versus the billiard shape parameter 
A for the Robnik billiard with a fixed density of boundary discretization 6 = 12. The numerical 
points are denoted by •, which are joined by a line just to guide the eye. 

FIG. 3. Ensemble-averaged (over 100 lowest odd-odd eigenstates) absolute error (measured in 
units of the mean level spacing) (in logarithmic units) versus the billiard shape parameter a/R for 
the Bunimovich stadium with a fixed density of boundary discretization b = 12. The numerical 
points are denoted by •, which are joined by a line just to guide the eye. 

FIG. 4. Ensemble-averaged (over 100 lowest eigenstates) absolute error (measured in units 
of the mean level spacing) versus the billiard parameter R (the radius of inner circle) for the 
desymmetrized Sinai billiard with a fixed density of boundary discretization b = 12. The numerical 
points are denoted by •, which are joined by a line just to guide the eye. 

FIG. 5. Absolute error of eigenstates (in logarithmic units) versus eigenenergy for the lowest 
400 odd eigenstates of the Robnik billiard at A = 0.27. The averages over consecutive stretches (of 
100 states each) are given in TABLE II, demonstrating that (|A.E|) is quite stationary. 

FIG. 6. Notation of the angles and chords used in boundary integral method (BIM). 

FIG. 7. BIM error (measured in units of the mean level spacing) of eigenstates versus energy. 
The error is the difference between the BIM value and the exact value. About 95% of the lowest 
1000 odd states are shown. Plot (a) is for the full circle billiard and (b) for the cardioid billiard 
(A = 1/2). In both cases b is fixed, 6 = 6. 

FIG. 8. Ensemble-averaged (over about 95% of the lowest 100 odd eigenstates) absolute BIM 
error versus the density of boundary discretization b and the best power-law fit for the various 
billiards. • represents the numerical data and the curve is the best power law fit whose a and 
coefficient A are given in each box. In (a) we have a full circle billiard, in (b) the 3ir/2 segment of 
a circular annulus with inner and outer radia 0.45 and 0.5, respectively, in (c) and (d) the Robnik 
billiard with A = 0.15andl/2, respectively, in (e) the 1/4 2x2 Bunimovich stadium, and in (f) the 
1/4 Sinai billiard with radius 1/2 inside a square of size 2. 

FIG. 9. Power-law exponent a versus the billiard shape parameter A for the Robnik billiard. 
The error bars denote the standard deviation from the best fit. 
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FIG. 10. Ensemble-averaged absolute BIM error (in units of the mean level spacing) (\AE\) 
(in logarithmic units) with fixed b = 12 against the biliiard shape parameter A for the Robnik 
billiard. 
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